Proximal Distance Algorithms: Theory and Practice

Proximal distance algorithms combine the classical penalty method of constrained minimization with distance majorization. If f(x) is the loss function, and C is the constraint set in a constrained minimization problem, then the proximal distance principle mandates minimizing the penalized loss f(x)+ρ2dist(x,C)2 and following the solution xρ to its limit as ρ tends to ∞. At each iteration the squared Euclidean distance dist(x,C)2 is majorized by the spherical quadratic ‖x− PC(xk)‖2, where PC(xk) denotes the projection of the current iterate xk onto C. The minimum of the surrogate function f(x)+ρ2‖x−PC(xk)‖2 is given by the proximal map proxρ−1f[PC(xk)]. The next iterate xk+1 automatically decreases the original penalized loss for fixed ρ. Since many explicit projections and proximal maps are known, it is straightforward to derive and implement novel optimization algorithms in this setting. These algorithms can take hundreds if not thousands of iterations to converge, but the simple nature of each iteration makes proximal distance algorithms competitive with traditional algorithms. For convex problems, proximal distance algorithms reduce to proximal gradient algorithms and therefore enjoy well understood convergence properties. For nonconvex problems, one can attack convergence by invoking Zangwill’s theorem. Our numerical examples demonstrate the utility of proximal distance algorithms in various high-dimensional settings, including a) linear programming, b) constrained least squares, c) projection to the closest kinship matrix, d) projection onto a second-order cone constraint, e) calculation of Horn’s copositive matrix index, f) linear complementarity programming, and g) sparse principal components analysis. The proximal distance algorithm in each case is competitive or superior in speed to traditional methods such as the interior point method and the alternating direction method of multipliers (ADMM). Source code for the numerical examples can be found at https://github.com/klkeys/proxdist.


Introduction
The solution of constrained optimization problems is part science and part art. As mathematical scientists explore the largely uncharted territory of high-dimensional nonconvex problems, it is imperative to consider new methods. The current paper studies a class of optimization algorithms that combine Courant's penalty method of optimization (Beltrami, 1970;Courant, 1943) with the notion of a proximal operator (Bauschke and Combettes, 2011;Moreau, 1962;Parikh and Boyd, 2013). The classical penalty method turns constrained minimization of a function f(x) over a closed set C into unconstrained minimization. The general idea is to seek the minimum point of a penalized version f(x) +ρq(x) of f(x), where the penalty q(x) is nonnegative and vanishes precisely on C. If one follows the solution vector x ρ as ρ tends to ∞, then in the limit one recovers the constrained solution. The penalties of choice in the current paper are squared Euclidean distances dist(x,C) 2 = inf y∈C ‖x−y‖ 2 .
The formula defines the proximal map of a function f(x). Here ‖ · ‖ is again the standard Euclidean norm, and f(x) is typically assumed to be closed and convex. Projection onto a closed convex set C is realized by choosing f(x) to be the 0/∞ indicator δ C (x) of C. It is possible to drop the convexity assumption if f(x) is nonnegative or coercive. In so doing, prox f (y) may become multi-valued. For example, the minimum distance from a nonconvex set to an exterior point may be attained at multiple boundary points. The point x in the definition (1) can be restricted to a subset S of Euclidean space by replacing f(x) by f(x) + δ S (x), where δ S (x) is the indicator of S.
One of the virtues of exploiting proximal operators is that they have been thoroughly investigated. For a large number of functions f(x), the map prox cf (y) for c > 0 is either given by an exact formula or calculable by an efficient algorithm. The known formulas tend to be highly accurate. This is a plus because the classical penalty method suffers from ill conditioning for large values of the penalty constant. Although the penalty method seldom delivers exquisitely accurate solutions, moderate accuracy suffices for many problems.
including projected Landweber, alternating projection onto the intersection of two or more closed convex sets, the alternating-direction method of multipliers (ADMM), and fast iterative shrinkage thresholding algorithms (FISTA) (Beck and Teboulle, 2009;Combettes and Pesquet, 2011;Landweber, 1951). Applications of distance majorization are more recent (Chi et al., 2014;Lange and Keys, 2014;Xu et al., 2017). The overall strategy consists of replacing the distance penalty dist(x,C) 2 by the spherical quadratic ‖x − y k ‖ 2 , where y k is the projection of the kth iterate x k onto C. To form the next iterate, one then sets The MM (majorization-minimization) principle guarantees that x k+1 decreases the penalized loss. We call the combination of Courant's penalty method with distance majorization the proximal distance principle. Algorithms constructed according to the principle are proximal distance algorithms.
The current paper extends and deepens our previous preliminary treatments of the proximal distance principle. Details of implementation such as Nesterov acceleration matter in performance. We have found that squared distance penalties tend to work better than exact penalties. In the presence of convexity, it is now clear that every proximal distance algorithm reduces to a proximal gradient algorithm. Hence, convergence analysis can appeal to a venerable body of convex theory. This does not imply that the proximal distance algorithm is limited to convex problems. In fact, its most important applications may well be to nonconvex problems. A major focus of this paper is on practical exploration of the proximal distance algorithm.
In addition to reviewing the literature, the current paper presents some fresh ideas. Among the innovations are: a) recasting proximal distance algorithms with convex losses as concave-convex programs, b) providing new perspectives on convergence for both convex and nonconvex proximal distance algorithms, c) demonstrating the virtue of folding constraints into the domain of the loss, and d) treating in detail seven interesting examples. It is noteworthy that some our new convergence theory is pertinent to more general MM algorithms.
It is our sincere hope to enlist other mathematical scientists in expanding and clarifying this promising line of research. The reviewers of the current paper have correctly pointed out that we do not rigorously justify our choices of the penalty constant sequence ρ k . The recent paper by Li et al. (2017) may be a logical place to start in filling this theoretical gap. They deal with the problem of minimizing f(x) subject to Ax = b through the quadratic penalized objective f (x) + ρ 2 ‖Ax − b‖ 2 . For the right choices of the penalty sequence ρ k , their proximal gradient algorithm achieves a O(k −1 ) rate of convergence for f(x) strongly convex. As a substitute, we explore the classical problem of determining how accurately the solution y ρ of the problem min x f (x) + ρ 2 that vanishes precisely on C. Polyak's proof relies on strong differentiability assumptions. Our proof for the case q(x) = dist(x,C) relies on convexity and is much simpler.
As a preview, let us outline the remainder of our paper. Section 2 briefly sketches the underlying MM principle. We then show how to construct proximal distance algorithms from the MM principle and distance majorization. The section concludes with the derivation of a few broad categories proximal distance algorithms. Section 3 covers convergence theory for convex problems, while Section 4 provides a more general treatment of convergence for nonconvex problems. To avoid breaking the flow of our exposition, all proofs are relegated to the Appendix. Section 5 discusses our numerical experiments on various convex and nonconvex problems. Section 6 closes by indicating some future research directions.

Derivation
The derivation of our proximal distance algorithms exploits the majorization-minimization (MM) principle (Hunter and Lange, 2004;Lange, 2010). In minimizing a function f(x), the MM principle exploits a surrogate function g(x | x k ) that majorizes f(x) around the current iterate x k . Majorization mandates both domination g(x | x k ) ≥ f(x) for all feasible x and tangency g(x k | x k ) = f(x k ) at the anchor x k . If x k+1 minimizes g(x | x k ), then the descent property f(x k+1 ) ≤ f(x k ) follows from the string of inequalities and equalities Clever selection of the surrogate g(x | x k+1 ) can lead to a simple algorithm with an explicit update that requires little computation per iterate. The number of iterations until convergence of an MM algorithm depends on how tightly g(x | x k ) hugs f(x). Constraint satisfaction is built into any MM algorithm. If maximization of f(x) is desired, then the objective f(x) should dominate the surrogate g(x | x k ) subject to the tangency condition. The next iterate x k+1 is then chosen to maximize g(x | x k ). The minorization-maximization version of the MM principle guarantees the ascent property.
The constraint set C over which the loss f(x) is minimized can usually be expressed as an using a convex combination of the squared distances. The neutral choice α i = g ρ for an irrelevant constant c k . If we put y k = ∑ i = 1 m α i P C i x k , then by definition the minimum of the surrogate g ρ (x | x k ) occurs at the proximal point (2) We call this MM algorithm the proximal distance algorithm. The penalty q(x) is generally smooth because at any point x where the projection P C (x) is single valued (Borwein and Lewis, 2006;Lange, 2016). This is always true for convex sets and almost always true for nonconvex sets. For the moment, we will ignore the possibility that P C (x) is multi-valued.
For the special case of projection of an external point z onto the intersection C of the closed sets C i , one should take f (x) = 1 2 ‖z − x‖ 2 . The proximal distance iterates then obey the explicit formula Linear programming with arbitrary convex constraints is another example. Here the loss is f(x) = v t x, and the update reduces to If the proximal map is impossible to calculate, but f(x) is L-smooth (∇f(x) is Lipschitz with constant L), then one can substitute the standard majorization . Minimizing the sum of the loss majorization plus the penalty majorization leads to the MM update This is a gradient descent algorithm without an intervening proximal map.
In moderate-dimensional problems, local quadratic approximation of f(x) can lead to a viable algorithm. For instance, in generalized linear statistical models, Xu et al. (2017) suggest replacing the observed information matrix by the expected information matrix. The latter matrix has the advantage of being positive semidefinite. In our notation, if A k ≈ d 2 f(x k ), then an approximate quadratic surrogate is The natural impulse is to update x by the Newton step This choice does not necessarily decrease f(x).
Step halving or another form of backtracking restores the descent property.
A more valid concern is the effort expended in matrix inversion. If A k is dense and constant, then extracting the spectral decomposition VDV t of A reduces formula (4) to when ρ is sufficiently large. For a generalized linear model, parameter updating involves solving the linear system for W k a diagonal matrix with positive diagonal entries. This task is equivalent to minimizing the least squares criterion In the unweighted case, extracting the singular value decomposition Z = USV T facilitates solving the system of equations (5). The svd decomposition is especially cheap if there is a substantial mismatch between the number rows and columns of Z. For sparse Z, the conjugate gradient algorithm adapted to least squares (Paige and Saunders, 1982b) is subject to much less ill conditioning than the standard conjugate gradient algorithm. Indeed, the algorithm LSQR and its sparse version LSMR (Fong and Saunders, 2011) perform well even when the matrix Z t W k 1/2 , ρI t is ill conditioned.
The proximal distance principle also applies to unconstrained problems. For example, consider the problem of minimizing a penalized loss ℓ(x)+p(Ax). The presence of the linear transformation Ax in the penalty complicates optimization. The strategy of parameter splitting introduces a new variable y and minimizes ℓ(x) + p(y) subject to the constraint y = Ax. If P M (z) denotes projection onto the manifold then the constrained problem can be solved approximately by minimizing the function for large ρ. If P M (z k ) consists of two subvectors u k and v k corresponding to x k and y k , then the proximal distance updates are x k + 1 = prox ρ −1 u k and y k + 1 = prox ρ −1 p Given the matrix A is n × p, one can attack the projection by minimizing the function This leads to the solution If n < p, then the Woodbury formula reduces the expense of matrix inversion.
Traditionally, convex constraints have been posed as inequalities C = {x : a(x) ≤ t}. Parikh and Boyd (2013) point out how to project onto such sets. The relevant Lagrangian for projecting an external point y amounts to with λ ≥ 0. The corresponding stationarity condition can be interpreted as a[prox λa (y)] = t. One can solve this one-dimensional equation for λ by bisection. Once λ is available, x = prox λa (y) is available as well. Parikh and Boyd (2013) note that the value a[prox λa (y)] is decreasing in λ. One can verify their claim by implicit differentiation of equation (7). This gives and consequently the chain rule inequality d dλ a prox λa y = − da x I + λd

Convergence: Convex Case
In the presence of convexity, the proximal distance algorithm reduces to a proximal gradient algorithm. This follows from the representation involving the penalty q(x). Thus, the proximal distance algorithm can be expressed as In this regard, there is the implicit assumption that q(x) is 1-smooth. This is indeed the case.
According to the Moreau decomposition (Bauschke and Combettes, 2011), for a single closed convex set C Because proximal operators of closed convex functions are nonexpansive (Bauschke and Combettes, 2011), the result follows for a single set. For the general penalty q(x) with m sets, the Lipschitz constants are scaled by the convex coefficients α i and added to produce an overall Lipschitz constant of 1.
It is enlightening to view the proximal distance algorithm through the lens of concaveconvex programming. Recall that the function is closed and convex for any nonempty closed set C. Danskin's theorem (Lange, 2016) justifies the directional derivative expression This equality allows us to identify the subdifferential ∂s(x) as the convex hull convP C (x). For any y ∈ ∂s(x k ), the supporting hyperplane inequality entails where d is a constant not depending on x. The same majorization can be generated by rearranging the majorization when y is the convex combination ∑ i β i p i of vectors p i from P C (x k ). These facts demonstrate that the proximal distance algorithm minimizing ‖x‖ 2 is often strongly convex regardless of whether f(x) itself is convex. If we replace the penalty dist(x,C) 2 by the penalty dist(Dx,C) 2 for a matrix D, then the function s(Dx) is still closed and convex, and minimization of f (x) + ρ 2 dist(Dx, C) 2 can also be viewed as an exercise in concave-convex programming.
In the presence of convexity, the proximal distance algorithm is guaranteed to converge. Our exposition relies on well-known operator results (Bauschke and Combettes, 2011). Proximal operators in general and projection operators in particular are nonexpansive and averaged. By definition an averaged operator is a convex combination of a nonexpansive operator N(x) and the identity operator I. The averaged operators on ℝ p with α ∈ (0,1) form a convex set closed under functional composition. Furthermore, M(x) and the base operator N(x) share their fixed points. The celebrated theorem of Krasnosel'skii (1955) and Mann (1953) says that if an averaged operator M(x) = αx + (1 − α)N(x) possesses one or more fixed points, then the iteration scheme x k+1 = M(x k ) converges to a fixed point.
Given the choice y k = ∑ i = 1 for all x. In the presence of convexity, this is equivalent to the directional derivative inequality for all v, which is in turn equivalent to z minimizing h ρ (x). Hence, if h ρ (x) attains its minimum value, then the proximal distance iterates converge to a minimum point.
Convergence of the overall proximal distance algorithm is tied to the convergence of the classical penalty method (Beltrami, 1970). In our setting, the loss is f(x), and the penalty is Assuming the objective f(x) + ρq(x) is coercive for all ρ ≥ 0, the theory mandates that the solution path x ρ is bounded and any limit point of the path attains the minimum value of f(x) subject to the constraints. Furthermore, if f(x) is coercive and possesses a unique minimum point in the constraint set C, then the path x ρ converges to that point.
Proximal distance algorithms often converge at a painfully slow rate. Following Mairal (2013), one can readily exhibit a precise bound.

Proposition 1 Suppose C is closed and convex and f(x) is convex. If the point z minimizes
suggests that Nesterov acceleration may vastly improve the chances for convergence. Nesterov acceleration for the general proximal gradient algorithm with loss ℓ(x) and penalty p(x) takes the form where L is the Lipschitz constant for ∇p(x) and d is typically chosen to be 3. Nesterov acceleration achieves an O(k −2 ) convergence rate (Su et al., 2014), which is vastly superior to the O(k −1 ) rate achieved by proximal gradient descent. The Nesterov update possesses the further desirable property of preserving affine constraints. In other words, if Ax k−1 = b and Ax k = b, then Az k = b as well.
In subsequent examples, we will accelerate our proximal distance algorithms by applying the algorithm map M(x) given by equation (2) to the shifted point z k of equation (10), yielding the accelerated update x k+1 = M(z k ). Algorithm 1 provides a schematic of a proximal distance algorithm with Nesterov acceleration. The recent paper of Ghadimi and Lan (2015) extends Nestorov acceleration to nonconvex settings.
In ideal circumstances, one can prove linear convergence of function values in the framework of Karimi et al. (2016).
Proposition 2 Suppose C is closed and convex and f(x) is L-smooth and μ-strongly convex.
Then h ρ (x) = f (x) + ρ 2 dist x, C 2 possesses a unique minimum point y, and the proximal We now turn to convergence of the penalty function iterates as the penalty constants ρ k tends to ∞. To simplify notation, we restrict attention to a single closed constraint set S. Let us start with a proposition requiring no convexity assumptions.
Proposition 3 If f(x) is continuous and coercive and S is compact, then the proximal distance iterates x k are bounded and the distance to the constraint set satisfies for some constant c. If in addition f(x) is continuously differentiable, then for some further constant d. Similar claims hold for the solutions y k of the penalty problem x, S 2 except that the assumption that S is compact can be dropped.
As a corollary, if the penalty sequence ρ k tends to ∞, then all limit points of x k must obey the constraint. Proposition 3 puts us into position to prove the next important result.
Proposition 4 If f(x) is continuously differentiable and coercive and S is convex, then the penalty function iterates defined by where y attains the constrained minimum and d is the constant identified in Proposition 3.

Convergence: General Case
Our strategy for addressing convergence in nonconvex problems fixes ρ and relies on Zangwill's global convergence theorem (Luenberger and Ye, 1984). This result depends in turn on the notion of a closed multi-valued map N(x). If x k converges to x ∞ and y k ∈ N(x k ) converges to y ∞ , then for N(x) to be closed, we must have y ∞ ∈ N(x ∞ ). The next proposition furnishes a prominent example.
Zangwill's global convergence theorem is phrased in terms of an algorithm map M(x) and a real-valued objective h(x). The theorem requires a critical set Γ outside which M(x) is closed. Furthermore, all iterates x k+1 ∈ M(x k ) must fall within a compact set. Finally, the descent condition h(y) ≤ h(x) should hold for all y ∈ M(x), with strict inequality when x ∉ Γ.
If these conditions are valid, then every convergent subsequence of x k tends to a point in Γ.
In the proximal distance context, we define the complement of Γ to consist of the points x with for all y ∈ M(x). This definition plus the monotonic nature of the proximal distance In general, the algorithm map M(x) is multi-valued in two senses. First, for a given z k ∈ P S (x k ), the minimum may be achieved at multiple points. This contingency is ruled out if the proximal map of f(x) is unique. Second, because S may be nonconvex, the projection may be rare. Accordingly, it makes no practical difference that we restrict the anchor points z k to lie in P S (x k ) rather than in convP S (x k ).
Proposition 6 If S is a closed nonempty set in ℝ p , then the projection operator P S (x) is single valued except on a set of Lebesgue measure 0.
In view of the preceding results, one can easily verify the next proposition.

Proposition 7
The algorithm map M(x) is everywhere closed.
To apply Zangwill's global convergence theory, we must in addition prove that the iterates x k+1 = M(x k ) remain within a compact set. This is true whenever the objective is coercive since the algorithm is a descent algorithm. As noted earlier, the coercivity of f(x) is a sufficient condition. One can readily concoct other sufficient conditions. For example, if f(x) is bounded below, say nonnegative, and S is compact, then the objective is also coercive. Indeed, if S is contained in the ball of radius r about the origin, then which proves that dist(x,S) is coercive. The next proposition summarizes these findings.
Proposition 8 If S is closed and nonempty, the objective f (x) + 1 2 dist x, S 2 is coercive, and the proximal operator prox ρ − 1f (x) is everywhere nonempty, then all limit points of the iterates x k+1 ∈ M(x k ) of the proximal distance algorithm occur in the critical set Γ.
This result is slightly disappointing. A limit point x could potentially exist with improvement in the objective for some but not all y ∈ convP S (x). This fault is mitigated by the fact that P S (x) is almost always single valued. In common with other algorithms in nonconvex optimization, we also cannot rule out convergence to a local minimum or a saddlepoint. One can improve on Proposition 8 by assuming that the surrogates g ρ (x | x k ) are all μ-strongly convex. This is a small concession to make because ρ is typically large. If f(x) is convex, then g ρ (x | x k ) is ρ-strongly convex by definition. It is also worth noting that any convex MM surrogate g(x | x k ) can be made μ-strongly convex by adding the viscosity penalty μ 2 ‖x − x k ‖ 2 majorizing 0. The addition of a viscosity penalty seldom complicates finding the next iterate x n+1 and has little impact on the rate of convergence when μ > 0 is small.

Proposition 9
Under the μ-strongly convexity assumption on the surrogates g ρ (x | x k ), the proximal distance iterates satisfy lim k→∞ ‖x k+1 −x k ‖ = 0. As a consequence, the set of limit points is connected as well as closed. Furthermore, if each limit point is isolated, then the iterates converge to a critical point.
2003). If h(x) is a function mapping ℝ p into ℝ ∪ +∞ , then its Fréchet subdifferential at x ∈ dom f is the set The set ∂ F h(x) is closed, convex, and possibly empty. If h(x) is convex, then ∂ F h(x) reduces to its convex subdifferential. If h(x) is differentiable, then ∂ F h(x) reduces to its ordinary differential. At a local minimum x, Fermat's rule 0 ∈ ∂ F h(x) holds.

Proposition 10
In an MM algorithm, suppose that h(x) is coercive, that the surrogates g(x | x k ) are differentiable, and that the algorithm map M(x) is closed. Then every limit point z of the MM sequence x k is critical in the sense that 0 ∈ ∂ F (−h)(z).
We will also need to invoke Łojasiewicz's inequality. This deep result depends on some rather arcane algebraic geometry (Bierstone and Milman, 1988;Bochnak et al., 2013). It applies to semialgebraic functions and their more inclusive cousins semianalytic functions and subanalytic functions. For simplicity we focus on semialgebraic functions. The class of semialgebraic subsets of ℝ p is the smallest class that: a. contains all sets of the form {x : q(x) > 0} for a polynomial q(x) in p variables,

b.
is closed under the formation of finite unions, finite intersections, and set complementation.
A function a : ℝ p ℝ r is said to be semialgebraic if its graph is a semialgebraic set of ℝ p + r . The class of real-valued semialgebraic functions contains all polynomials p(x) and all 0/1 indicators of algebraic sets. It is closed under the formation of sums, products, absolute values, reciprocals when a(x) 6≠ 0, nth roots when a(x) ≥ 0, and maxima max{a(x), b(x)} and minima min{a(x), b(x)}. For our purposes, it is important to note that dist(x,S) is a semialgebraic function whenever S is a semialgebraic set.
Łojasiewicz's inequality in its modern form (Bolte et al., 2007) requires a function h(x) to be closed (lower semicontinuous) and subanalytic with a closed domain. If z is a critical point of h(x), then Here the exponent θ ∈ [0,1), the radius r, and the constant c depend on z. This inequality is valid for semialgebraic functions since they are automatically subanalytic. We will apply Łojasiewicz's inequality to the limit points of an MM algorithm. The next proposition is an elaboration and expansion

Proposition 11
In an MM algorithm suppose the objective h(x) is coercive, continuous, and subanalytic and all surrogates g(x | x k ) are continuous, μ-strongly convex, and satisfy the Lsmoothness condition converge to a critical point.
The last proposition applies to proximal distance algorithms. The loss f(x) must be subanalytic and differentiable with a locally Lipschitz gradient. Furthermore, all surrogates ‖x − y k ‖ 2 should be coercive and μ-strongly convex. Finally, the constraints sets S i should be subanalytic. Semialgebraic sets and functions will do. Under these conditions and regardless of how the projected points P Si (x) are chosen, the MM iterates are guaranteed to converge to a critical point.

Examples
The following examples highlight the versatility of proximal distance algorithms in a variety of convex and nonconvex settings. Programming details matter in solving these problems. Individual programs are not necessarily long, but care must be exercised in projecting onto constraints, choosing tuning schedules, folding constraints into the domain of the loss, implementing acceleration, and declaring convergence.
where ϵ 1 =10 −6 and ϵ 2 =10 −4 are typical values. The number of iterations until convergence is about 1000 in most examples. This handicap is offset by the simplicity of each stereotyped update. Our code is available as supplementary material to this paper. Readers are encouraged to try the code and adapt it to their own examples.

Linear Programming
Two different tactics suggest themselves for constructing a proximal distance algorithm. The first tactic rolls the standard affine constraints Ax = b into the domain of the loss function v t x. The standard nonnegativity requirement x ≥ 0 is achieved by penalization. Let x k be the current iterate and y k = (x k ) + be its projection onto ℝ + n . Derivation of the proximal distance algorithm relies on the Lagrangian v t x + ρ 2 ‖x − y k ‖ 2 + λ One can multiply the corresponding stationarity equation by A and solve for the Lagrange multiplier λ in the form assuming A has full row rank. Inserting this value into the stationarity equation gives the MM update where A − = A t (AA t ) −1 is the pseudo-inverse of A.
The second tactic folds the nonnegativity constraints into the domain of the loss. Let p k denote the projection of x k onto the affine constraint set Ax = b. Fortunately, the surrogate ‖x − p k ‖ 2 splits the parameters. Minimizing one component at a time gives the update x k+1 with components The projection p k can be computed via where A − is again the pseudo-inverse of A. , which relies on a fast implementation of ADMM. The second is the commercial Gurobi solver, which ships with implementations of both the simplex method and a barrier (interior point) method; in this example, we use its barrier algorithm. The first seven rows of the table summarize linear programs with dense data A, b, and v. The bottom six rows rely on random sparse matrices A with sparsity level 0.01. For dense problems, the proximal distance algorithms start the penalty constant ρ at 1 and double it every 100 iterations.
Because we precompute and cache the pseudoinverse A − of A, the updates (12) and (13) reduce to vector additions and matrix-vector multiplications.
For sparse problems the proximal distance algorithms update ρ by a factor of 1.5 every 50 iterations. To avoid computing large pseudoinverses, we appeal to the LSQR variant of the conjugate gradient method (Paige and Saunders, 1982b,a) to solve the linear systems (11) and (14). The optima of all four methods agree to about 4 digits of accuracy. It is hard to declare an absolute winner in these comparisons. Gurobi and SCS clearly perform better on low-dimensional problems, but the proximal distance algorithms are competitive as dimensions increase. PD1, the proximal distance algorithm over an affine domain, tends to be more accurate than PD2. If high accuracy is not a concern, then the proximal distance algorithms are easily accelerated with a more aggressive update schedule for ρ.

Constrained Least Squares
Constrained least squares programming subsumes constrained quadratic programming. A typical quadratic program involves minimizing the quadratic 1 2 x t Qx − p t x subject to x ∈ C for a positive definite matrix Q. Quadratic programming can be reformulated as least squares by taking the Cholesky decomposition Q = LL t of Q and noting that The constraint x ∈ C applies in both settings. It is particularly advantageous to reframe a quadratic program as a least squares problem when Q is already presented in factored form or when it is nearly singular (Bemporad, 2018). To simplify subsequent notation, we replace L t by the rectangular matrix A and L −1 p by y. The key to solving constrained least squares is to express the proximal distance surrogate as as in equation (6). As noted earlier, in sparse problems the update x k+1 can be found by a fast stable conjugate gradient solver. Table 2 compares the performance of the proximal distance algorithm for least squares estimation with probability-simplex constraints to the open source nonlinear interior point Simplex constrained problems arise in hyperspectral imaging (Heylen et al., 2011;Keshava, 2003), portfolio optimization (Markowitz, 1952), and density estimation (Bunea et al., 2010). Test problems were generated by filling an n×p matrix A and an n-vector y with standard normal deviates. For sparse problems we set the sparsity level of A to be 10/p. Our setup ensures that A has full rank and that the quadratic program has a solution. For the proximal distance algorithm, we start ρ at 1 and multiply it by 1.5 every 200 iterations. Table   2 suggests that the proximal distance algorithm and the interior point solvers perform equally well on small dense problems. However, in high-dimensional and low-accuracy environments, the proximal distance algorithm provides much better scalability.

Closest Kinship Matrix
In genetics studies, kinship is measured by the fraction of genes two individuals share identical by descent. For a given pedigree, the kinship coefficients for all pairs of individuals appear as entries in a symmetric kinship matrix Y. This matrix possesses three crucial properties: a) it is positive semidefinite, b) its entries are nonnegative, and c) its diagonal entries are 1 2 unless some pedigree members are inbred. Inbreeding is the exception rather than the rule. Kinship matrices can be estimated empirically from single nucleotide polymorphism (SNP) data, but there is no guarantee that the three highlighted properties are satisfied. Hence, it helpful to project Y to the nearest qualifying matrix.
This projection problem is best solved by folding the positive semidefinite constraint into the domain of the Frobenius loss function 1 2 ‖X − Y‖ F 2 . As we shall see, the alternative of imposing two penalties rather than one is slower and less accurate. Projection onto the constraints implied by conditions b) and c) is trivial. All diagonal entries x ii of X are reset to 1 2 , and all off-diagonal entries x ij are reset to max{x ij ,0}. If P(X k ) denotes the current projection, then the proximal distance algorithm minimizes the surrogate where c k is an irrelevant constant. The minimum is found by extracting the spectral decomposition UDU t of 1 1 + ρ Y + ρ 1 + ρ P X k and truncating the negative eigenvalues. This gives the update X k+1 = UD + U t in obvious notation. This proximal distance algorithm and its Nesterov acceleration are simple to implement in a numerically oriented language such as Julia. The most onerous part of the calculation is clearly the repeated eigen-decompositions. Table 3 compares three versions of the proximal distance algorithm to Dykstra's algorithm (Boyle and Dykstra, 1986). Higham proposed Dykstra's algorithm for the related problem of finding the closest correlation matrix Higham (2002). In Table 3 algorithm PD1 is the unadorned proximal distance algorithm, PD2 is the accelerated proximal distance, and PD3 is the accelerated proximal distance algorithm with the positive semidefinite constraints folded into the domain of the loss. On this demanding problem, these algorithms are comparable to Dykstra's algorithm in speed but slightly less accurate. Acceleration of the proximal distance algorithm is effective in reducing both execution time and error. Folding the positive semidefinite constraint into the domain of the loss function leads to further improvements. The data matrices M in these trials were populated by standard normal deviates and then symmetrized by averaging opposing off-diagonal entries. In algorithm PD1 we set ρ k = max{1.2 k ,2 22 }. In the accelerated versions PD2 and PD3 we started ρ at 1 and multiplied it by 5 every 100 iterations. At the expense of longer compute times, better accuracy can be achieved by all three proximal distance algorithms with a less aggressive update schedule.

Projection onto a Second-Order Cone Constraint
Second-order cone programming is one of the unifying themes of convex analysis (Alizadeh and Goldfarb, 2003;Lobo et al., 1998). It revolves around conic constraints of the form {u : ‖Au + b‖ ≤ c t u + d}. Projection of a vector x onto such a constraint is facilitated by parameter splitting. In this setting parameter splitting introduces a vector w, a scalar r, and the two affine constraints w = Au + b and r = c t u + d. The conic constraint then reduces to the Lorentz cone constraint ‖w‖ ≤ r, for which projection is straightforward (Boyd and Taking a cue from Example 5.1, we incorporate the affine constraints in the domain of the objective function. If we represent projection onto L by then the Lagrangian generated by the proximal distance algorithm amounts to Solving for the multipliers λ and θ in equations (16) and (17) and substituting their values in equation (15) yield This leads to the MM update The updates w k+1 = Au k+1 + b and r k+1 = c t u k+1 + d follow from the constraints. Table 4 compares the proximal distance algorithm to SCS and Gurobi. Echoing previous examples, we tailor the update schedule for ρ differently for dense and sparse problems. Dense problems converge quickly and accurately when we set ρ 0 = 1 and double ρ every 100 iterations. Sparse problems require a greater range and faster updates of ρ, so we set ρ 0 = 0.01 and then multiply ρ by 2.5 every 10 iterations. For dense problems, it is clearly With a large sparse constraint matrix A, extraction of its spectral decomposition becomes prohibitive. If we let E = (ρ −1/2 I A t c), then we must solve a linear system of equations defined by the Gramian matrix G = EE t . There are three reasonable options for solving this system. The first relies on computing and caching a sparse Cholesky decomposition of G. The second computes the QR decomposition of the sparse matrix E. The R part of the QR decomposition coincides with the Cholesky factor. Unfortunately, every time ρ changes, the Cholesky or QR decomposition must be redone. The third option is the conjugate gradient algorithm. In our experience the QR decomposition offers superior stability and accuracy. When E is very sparse, the QR decomposition is often much faster than the Cholesky decomposition because it avoids forming the dense matrix A t A. Even when only 5% of the entries of A are nonzero, 90% of the entries of A t A can be nonzero. If exquisite accuracy is not a concern, then the conjugate gradient method provides the fastest update. Table 4 reflects this choice. The value μ(M) = 0 is attained for the vectors 1 2 odd dimension cannot be written as a positive semidefinite matrix, a nonnegative matrix, or a sum of two such matrices.

Copositive Matrices
The proximal distance algorithm minimizes the criterion and generates the updates It takes a gentle tuning schedule to get decent results. The choice ρ k = 1.2 k converges in 600 to 700 iterations from random starting points and reliably yields objective values below 10 −5 for Horn matrices. The computational burden per iteration is significantly eased by exploiting the cached spectral decomposition of M. Table 5 compares the performance of the proximal distance algorithm to the Mosek solver on a range of Horn matrices. Mosek uses semidefinite programming to decide whether M can be decomposed into a sum of a positive semidefinite matrix and a nonnegative matrix. If not, Mosek declares the problem infeasible. Nesterov acceleration improves the final loss for the proximal distance algorithm, but it does not decrease overall computing time.
Testing for copositivity is challenging because neither the loss function nor the constraint set is convex. The proximal distance algorithm offers a fast screening device for checking whether a matrix is copositive. On random 1000×1000 symmetric matrices M, the method invariably returns a negative index in less than two seconds of computing time. Because the vast majority of symmetric matrices are not copositive, accurate estimation of the minimum is not required. Table 6 summarizes a few random trials with lower-dimensional symmetric matrices. In higher dimensions, Mosek becomes non-competitive, and Nesterov acceleration is of dubious value.

Linear Complementarity Problem
The linear complementarity problem (Murty and Yu, 1988) consists of finding vectors x and y with nonnegative components such that x t y = 0 and y = Ax + b for a given square matrix A and vector b. The natural loss function is 1 2 ‖y − Ax − b‖ 2 . To project a vector pair (u,v) onto the nonconvex constraint set, one considers each component pair (u i , v i ) in turn. If u i ≥ max{v i ,0}, then the nearest pair (x,y) has components (x i , y i ) = (u i ,0). If v i ≥ max{u i ,0}, then the nearest pair has components (x i , y i ) = (0, v i ). Otherwise, (x i , y i ) = (0,0). At each iteration the proximal distance algorithm minimizes the criterion Substituting the value of y from the second equation into the first equation leads to the updates The linear system (19) can be solved in low to moderate dimensions by computing and caching the spectral decomposition of A t A and in high dimensions by the conjugate gradient method. Table 7 compares the performance of the proximal gradient algorithm to the Gurobi solver on some randomly generated problems.

Sparse Principal Components Analysis
Let X be an n × p data matrix gathered on n cases and p predictors. Assume the columns of X are centered to have mean 0. Principal component analysis (PCA) (Hotelling, 1933;Pearson, 1901) operates on the sample covariance matrix S = 1 n X t X. Here we formulate a proximal distance algorithm for sparse PCA (SPCA), which has attracted substantial interest in the machine learning community (Berthet and Rigollet, 2013b,a;D'Aspremont et al., 2007;Johnstone and Lu, 2009;Journée et al., 2010;Witten et al., 2009;Zou et al., 2006). According to a result of Ky Fan (Fan, 1949), the first q principal components (PCs) u 1 ,…,u q can be extracted by maximizing the function tr(U t SU) subject to the matrix constraint U t U = I q , where u i is the ith column of the p×q matrix U. This constraint set is called a Stiefel manifold. One can impose sparsity by insisting that any given column u i have at most r nonzero entries. Alternatively, one can require the entire matrix U to have at most r nonzero entries. The latter choice permits sparsity to be distributed non-uniformly across columns.
Extraction of sparse PCs is difficult for three reasons. First, the Stiefel manifold M q and both sparsity sets are nonconvex. Second, the objective function is concave rather than convex. Third, there is no simple formula or algorithm for projecting onto the intersection of the two constraint sets. Fortunately, it is straightforward to project onto each separately. Let P M q calculated by extracting a partial singular value decomposition U = VΣW t of U and setting P Mq (U) = VW t (Golub and Van Loan, 2012). Here V and W are orthogonal matrices of dimension p×q and q×q, respectively, and Σ is a diagonal matrix of dimension q × q. Let P S r (U) denote the projection of U onto the sparsity set S r = V : v i j ≠ 0 for at most r entries of each column v i .
Because P S r (U) operates column by column, it suffices to project each column vector u i to sparsity. This entails nothing more than sorting the entries of u i by magnitude, saving the r largest, and sending the remaining p−r entries to 0. If the entire matrix U must have at most r nonzero entries, then U can be treated as a concatenated vector during projection.
The key to a good algorithm is to incorporate the Stiefel constraints into the domain of the objective function (Kiers, 1990;Kiers and ten Berge, 1992) and the sparsity constraints into the distance penalty. Thus, we propose decreasing the criterion at each iteration subject to the Stiefel constraints. The loss can be majorized via  Shen and Huang, 2008). When the loading vectors of U are orthogonal, this criterion reduces to the familiar definition tr(U t X t XU)/tr(X t X) of PVE for ordinary PCA. The proximal distance algorithm enforces either matrix-wise or column-wise sparsity. In contrast, SPC enforces only column-wise sparsity via the constraint ‖u i ‖ 1 ≤ c for each column u i of U. We take c = 8. The number of nonzeroes per loading vector output by SPC dictates the sparsity level for the column-wise version of the proximal distance algorithm. Summing these counts across all columns dictates the sparsity level for the matrix version of the proximal distance algorithm.
Figures 1 and 2 demonstrate the superior PVE and computational speed of both proximal distance algorithms versus SPC. The type of projection does not appear to affect the computational performance of the proximal distance algorithm, as both versions scale equally well with q. However, the matrix projection, which permits the algorithm to more freely assign nonzeroes to the loadings, attains better PVE than the more restrictive columnwise projection. For both variants of the proximal distance algorithm, Nesterov acceleration improves both fitting accuracy and computational speed, especially as the number of PCs q increases.

Discussion
The proximal distance algorithm applies to a host of problems. In addition to the linear and quadratic programming examples considered here, our previous paper (Lange and Keys, 2014) derives and tests algorithms for binary piecewise-linear programming, ℓ 0 regression, matrix completion (Cai et al., 2010;Candès and Tao, 2010;Chen et al., 2012;Mazumder et al., 2010), and sparse precision matrix estimation (Friedman et al., 2008). Other potential applications immediately come to mind. An integer linear program in standard form can be expressed as minimizing c t x subject to Ax + s = b, s ≥ 0, and x ∈ ℤ p . The latter two constraints can be combined in a single constraint for which projection is trivial. The affine constraints should be folded into the domain of the objective. Integer programming is NP hard, so that the proximal distance algorithm just sketched is merely heuristic. Integer linear programming includes traditional NP hard problems such as the traveling salesman problem, the vertex cover problem, set packing, and Boolean satisfiability. It will be interesting to see if the proximal distance principle is competitive in meeting these challenges. Our experience with the closest lattice point problem (Agrell et al., 2002) and the eight queens problem suggests that the proximal distance algorithm can be too greedy for hard combinatorial problems. The nonconvex problems solved in this paper are in some vague sense easy combinatorial problems.
constraints. Incrementing ρ too quickly causes the algorithm to veer off the solution path guaranteed by the penalty method. Given the chance of roundoff error even with double precision arithmetic, it is unwise to take ρ all the way to ∞. Trial and error can help in deciding whether a given class of problems will benefit from an aggressive update schedule and strict or loose convergence criteria. In problems with little curvature such as linear programming, more conservative updates are probably prudent. The linear programming, closest kinship matrix, and SPCA problems document the value of folding constraints into the domain of the loss. In the same spirit it is wise to minimize the number of constraints. A single penalty for projecting onto the intersection of two constraint sets is almost always preferable to the sum of two penalties for their separate projections. Exceptions to this rule obviously occur when projection onto the intersection is intractable. The integer linear programming problem mentioned previously illustrates these ideas.
Our earlier proximal distance algorithms ignored acceleration. In many cases the solutions produced had very low accuracy. The realization that convex proximal distance algorithms can be phrased as proximal gradient algorithms convinced us to try Nesterov acceleration. We now do this routinely on the subproblems with ρ fixed. This typically forces tighter path following and a reduction in overall computing times. Our examples generally bear out the contention that Nesterov acceleration is useful in nonconvex problems (Ghadimi and Lan, 2015). It is noteworthy that the value of acceleration often lies in improving the quality of a solution as much as it does in increasing the rate of convergence. Of course, acceleration cannot prevent convergence to an inferior local minimum.
On both convex and nonconvex problems, proximal distance algorithms enjoy global convergence guarantees. On nonconvex problems, one must confine attention to subanalytic sets and subanalytic functions. This minor restriction is not a handicap in practice. Determining local convergence rates is a more vexing issue. For convex problems, we review existing theory for a fixed penalty constant ρ. We hope readers will sense the potential of the proximal distance principle. This simple idea offers insight into many existing algorithms and a straightforward path in devising new ones. Effective proximal and projection operators usually spell the difference between success and failure. The number and variety of such operators is expanding quickly as the field of optimization relinquishes it fixation on convexity. The current paper research leaves many open questions about tuning schedules, rates of convergence, and acceleration in the face of nonconvexity. We welcome the contributions of other mathematical scientists in unraveling these mysteries and in inventing new proximal distance algorithms.
h ρ Adding the result over k and invoking the descent property h ρ (x k+1 ) ≤ h ρ (x k ) produce the error bound Setting x equal to a minimal point z gives the stated result. ■

A.2. Proposition 2
The existence and uniqueness of y are obvious. The remainder of the proof hinges on the assumptions that h ρ (x) is μ-strongly convex and the surrogate g ρ (x | x k ) is L + ρ smooth. The latter assumption yields h ρ x − h ρ y ≤ g ρ x | y − g ρ y | y ≤ ∇g ρ y | y t x − y + L + ρ 2 ‖x − y‖ 2 = L + ρ 2 ‖x − y‖ 2 .
The strong convexity condition h ρ (y) − h ρ (x) ≥ ∇h ρ x t y − x + is valid by definition, where y attains the constrained minimum. Thus, coerciveness implies that the sequence y k is bounded. When f(x) is continuously differentiable, the proof of the second claim also applies if we substitute y k for x k and P S (y k ) for P S (x k−1 ). ■

A.4. Proposition 4
Because the function f (x) + ρ k 2 dist x, S 2 is convex and has the value f(y) and gradient ∇f(y) at a constrained minimum y, the supporting hyperplane principle says The first-order optimality condition ∇f(y) t [P S (y k ) − y] ≥ 0 holds given y is a constrained minimum. Hence, the Cauchy-Schwarz inequality and Proposition 3 imply

A.5. Proposition 5
Let x k converge to x ∞ and y k ∈ P S (x k ) converge to y ∞ . For an arbitrary y ∈ S, taking limits in the inequality ‖x k −y k ‖ ≤ ‖x k -y‖ yields ‖x ∞ −y ∞ ‖ ≤ ‖x ∞ −y‖; consequently, y ∞ ∈ P S (x ∞ ). To prove the second assertion, take y k ∈ P S (x k ) and observe that Rearranging this inequality and summing over k yield ∑ n = 0 Thus, the sequence x k is a fast Cauchy sequence and converges to a unique limit in W. ■ Proportion of variance explained by q PCs for each algorithm. Here PD1 is the accelerated proximal distance algorithm enforcing matrix sparsity, PD2 is the accelerated proximal distance algorithm enforcing column-wise sparsity, and SPC is the orthogonal sparse PCA method from PMA.  Computation times for q PCs for each algorithm. Here PD1 is the accelerated proximal distance algorithm enforcing matrix sparsity, PD2 is the accelerated proximal distance algorithm enforcing column-wise sparsity, and SPC is the orthogonal sparse PCA method from PMA.  CPU times and optima for the closest kinship matrix problem. Here the kinship matrix is n × n, PD1 is the proximal distance algorithm, PD2 is the accelerated proximal distance, PD3 is the accelerated proximal distance algorithm with the positive semidefinite constraints folded into the domain of the loss, and Dykstra is Dykstra's adaptation of alternating projections. All times are in seconds.